Cite this Article
Wang Bi-De, Song Jian, Li Ran, Han Ren, Zheng Gang, Yang Hui. A novel particle tracking velocimetry method for complex granular flow field. Chinese Physics B, 2020, 29(1): 014207  
Permissions
A novel particle tracking velocimetry method for complex granular flow field
Wang Bi-De1, Song Jian1, Li Ran2, Han Ren1, Zheng Gang2, Yang Hui1, †
School of Optical-Electrical and Computer Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China
School of Medical Instrument and Food Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China

 

† Corresponding author. E-mail: yanghui@usst.edu.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11572201 and 91634202).

Abstract

Particle tracking velocimetry (PTV) is one of the most commonly applied granular flow velocity measurement methods. However, traditional PTV methods may have issues such as high mismatching rates and a narrow measurement range when measuring granular flows with large bulk density and high-speed contrast. In this study, a novel PTV method is introduced to solve these problems using an optical flow matching algorithm with two further processing steps. The first step involves displacement correction, which is used to solve the mismatching problem in the case of high stacking density. The other step is trajectory splicing, which is used to solve the problem of a measurement range reduction in the case of high-speed contrast The hopper flow experimental results demonstrate superior performance of this proposed method in controlling the number of mismatched particles and better measuring efficiency in comparison with the traditional PTV method.

PACS: 42.30.-d;47.11.-j;47.57.Gc
Keyword:particle tracking velocimetry;optical flow;displacement correction;trajectory splicing
1. Introduction

In the past few decades, many new optical non-contact measurement methods have emerged owing to developments in optical technology. Among them, spatial filtering velocimetry (SFV),[1,2] particle image velocimetry (PIV),[3–6] and particle tracking velocimetry (PTV)[7–10] are the three most commonly employed methods for granular flow fields. For a nearly uniform flow velocity distribution, results obtained by these three methods are highly similar. However, if granular flows are of high concentration, an inconvenient drawback is exposed, namely, PIV lags far behind PTV in terms of the ability to recognize particle displacement.[11] In the meantime, if the granular flow has a large velocity contrast, the SFV method may cause large measurement errors, since the measurement results are averaged. Arguably, PTV is the best method to measure complex granular flows.

The PTV technique mainly includes particle recognition and particle matching. Particle recognition is affected by hardware conditions, as the particle recognition accuracy is positively correlated with the camera imaging quality. After determining the hardware parameters, the particle-matching method will become a key factor affecting the accuracy of the PTV technique. Usually, the PTV matching algorithm needs to primarily set the matching rules of tracer particles to establish a particle-matching calculation model. The matching rules that appear in the traditional PTV method[12] include the following: particle motion displacement needs to be smaller than the distance between the particle and other particles, the motion of the tracer particles in a small space needs to be consistent, etc. When measuring complex granular flow scenes, such as hopper granular flows, the matching rules depicted above are not perfectly applicable, because the velocity field distribution of the hopper flow is not only numerically different, but directions of movement of the particles at different positions are also contrasting.[13] If these rules are enforced in the matching algorithm, it is likely to lead to a significant reduction in particle-matching accuracy. In 2011, Balevicius et al.[14] employed the discrete element method to study particle motion in different types of hoppers. In 2018, Ma et al.[15] used the same method to study particle clogging within hoppers. The latest research addressing hopper granular flow was conducted by Zhang et al.,[16] who studied the applicability of Beverloo’s scaling law within hoppers in 2019. Apparently, due to the lack of measurement methods, many scholars prefer simulation to experiment.

The optical flow method is a motion detection method in the field of computer vision. Compared with other methods in this field, such as the CamShift algorithm,[17] the optical flow algorithm exhibits better performance in tracking multiple moving targets. Optical flow methods are usually divided into sparse optical flow and dense optical flow. The representative method of dense optical flows is the Horn–Schunck (H–S) method,[18] whereas the representative methods of sparse optical flows are the Lucas–Kanade (L–K) method[19] and Pyramid Lucas–Kanade (PyrLK) method.[20] Sparse optical flow is employed for image registration specifically for sparse points on an image, which means that this method needs to track certain coordinates within an image. Due to the small amount of calculation and high precision,[21] it is often used for image registration and target tracking. Dense optical flow requires the calculation of the offset of all the points in the image to form a dense field.[22] The advantage of this method is that the target can be registered at the pixel level, whereas the disadvantage is that it requires significantly more calculation than the sparse optical flow. In the field of granular flow velocity measurement, we need to measure the motion of a single particle, such that the sparse optical flow method is more appropriate than the dense optical flow method. Considering the granular flow movement in some complex scenarios, the PyrLK method exhibits a better performance in tracking particles with large displacements compared with the L–K method.

The PyrLK method uses corner displacement to characterize particles. However, in complex scenarios such as hopper granular flow, the centroid displacement provides a more accurate representation of particle motion. In this study, a displacement correction method of corner-to-centroid is proposed for the improvement of the PyrLK method. Theoretically, the trajectory reliability of the particle between two frames has a positive correlation with the camera frame rate. Therefore, to measure the particle motion over a longer time scale, it is necessary to splice all trajectories at this scale.

Based on the claim above, choosing hopper granular flow as a typical complex scenario, this study introduces an improved PTV method compounding displacement correction and trajectory splicing for the measurement of complex granular flows, and builds an experimental device. The results show that the designed method effectively improves both precision and efficiency compared with the traditional PTV method.

2. PyrLK-based PTV method

The basic principle of the PTV technique is to track the motion trajectory of the same particle at different time points, calculate the motion displacement of the particle using the trajectory, and finally obtain the motion velocity of the particle by the displacement and time interval.[23,24] Determining the position of the particle between two frames is a critical process in the PTV method. The PyrLK method determines the pixel displacement in the image sequence by measuring the change of the pixel intensity data in the image sequence and the correlation between adjacent frames,[18] thereby performing tracking of the particle target. Applying the PyrLK method to the PTV matching algorithm indicates that there is no need to assume matching rules for moving targets, which is highly advantageous for the measurement of complex flow fields.

The PTV matching algorithm employing the PyrLK method includes three steps: extracting corner points, establishing constraint equations, and layering by the Gaussian pyramid. The specific implementation method will be described in detail later in this paper.

2.1. Extracting corner points

The PyrLK method selects corner points as the particle feature points. If a small change in a point on an image causes a large change in the gray-scale, this point will be considered as a corner coordinate. A common method for corner tracking is the Shi–Tomashi feature detection, and the tracking principle is as follows.

Assume (x, y) as the coordinates of a pixel in an image, I(x, y) is the intensity information of this pixel. After the displacement of (Δx, Δy), the gray-scale variation function, E(Δx, Δy), of pixel (x, y) can be expressed as follows:

where w(x, y) represents the weighting factor of the pixel (x, y). Expand the function I (x + Δx, y + Δy) using the Taylor formula
Taking the first-order approximation expression of Eq. (2) and substituting it into Eq. (1), we can obtain
where matrix M represents a structural tensor, and its two eigenvalues λ1 and λ2 represent the change rates of I(x, y) in two directions. As long as λ1 and λ2 of a certain point are large, image movement in any direction will result in significant gray-scale changes, and we refer to this point as the corner point.

After the Shi–Tomashi feature detection, the obtained corner coordinate is an integer value due to the limitations of the acquisition camera, despite the fact that the corner coordinate position is rarely at the integer position. To overcome this problem, we decide to pursue sub-pixel corner detection techniques,[25,26] plot and fit curves based on the image pixel values, and use mathematical operations to determine where the peak appears between the pixels. Figure 1 is the result of the Shi–Tomashi feature detection employed to extract the corner coordinate in a granular flow image.

Fig. 1. (a) Original image of granular flow and (b) image using Shi–Tomashi feature detection to extract corner coordinates.
2.2. Establishing constraint equations

Suppose that the pixel coordinates (x, y) of the corner point have been determined, and the intensity information function in time t is I(x, y, t). After Δt time interval, the intensity information function is expressed as I(x + Δx, y + Δy, t + Δt) If the movement of this point is small and continuous[27] we can use the Taylor formula to expand I(x + Δx, y + Δy, t + Δt) as follows:

where ε is the second-order infinitesimal of Taylor expansion. When Δt approaches zero, the second-order infinitesimal of Taylor expansion can be ignored, and the intensity is assumed to be constant. At this time, the basic optical flow constraint equation can be expressed as follows:
where u = dx/dt represents the horizontal velocity component of the corner point, while w = dy/dt represents the vertical velocity component. As a single equation cannot solve a problem with two unknowns u and w, an assumption is made within the L–K method, that the optical flow in the adjacent region is consistent, which indicates that adjacent pixels on the same surface have similar motion. Through the assumption of spatial consistency, n optical flow constraint equations can be established using n pixels in the adjacent region to estimate the optical flow component An×2⋅d2×1=[ I1xI2x⋮InxI1yI2y⋮Iny ][ uw ]=−[ I1tI2t⋮Int ]=bn×1 where A is a two-dimensional matrix obtained by merging (∂I/∂x, ∂I/∂y) of n pixels, d is a one-dimensional matrix obtained by merging u and w, and b is a one-dimensional matrix obtained by merging (∂I/∂t) of n pixels.

The least squares method is used to find optimal solutions of these n equations, and the motion components of u and w can be obtained. The detailed calculation steps are shown by the following equations: (A TA)⋅d=[ ∑​IxIx∑​IxIy∑​IyIx∑​IyIy ][ uw ]=−[ ∑​IxIt∑​IyIt ]=(ATb),

2.3. Layering by Gaussian pyramid

The condition for the Taylor expansion in the second step is that the corner position does not change drastically. Once the moving target velocity increases, the displacement becomes larger, and the algorithm produces a larger error. To solve these drawbacks, this study proposes to use the PyrLK method, which is an improved L–K method based on the pyramid stratification proposed by Bouguet.[28] The principle behind the PyrLK method is explained in Fig. 2. First, the optical flow is calculated from the top layer of the image pyramid, and subsequently the calculation result of the top layer is used as the starting value of the calculation of the next layer, and the process is repeated until the bottom layer of the pyramid. After layer-by-layer processing, the movement of the corner point is small, and continuity is satisfied. Eventually, scenes with high speed and large displacements can be detected.

Fig. 2. Diagram of Gaussian pyramid layering.
3. Further processing
3.1. Displacement correction

In some scenarios, since velocities in the flow field are not completely the same, particles may also rotate during motion. In this case, centroid displacement is more representative of the actual motion than corner displacement, hence the corner coordinates need to be converted to centroid coordinates before the velocity calculation.

The specific displacement correction method used in this study is as follows. First, a reasonable threshold is determined by observing the segmented image until all tracer particles are highlighted. The image is then subjected to binarization processing, resulting in a processed image that contains the physical contour of the particles. After extracting contours of respective entities, a key step is to evaluate the affiliation between the corner point and contour of each particle. Finally, a geometric moment algorithm is used to obtain the centroid of the particle contour, and the centroid matching of the tracer particles is realized.

Evaluating the affiliation between the corner point and contour is key to the whole displacement correction process. The implementation of this step is explained in Fig. 3: for each corner point, a loop is set to determine the positional relationship between the contours of the entities in the image and the corner points. When the corner point is within the contour or on the contour edge, it can be directly determined that the corner point belongs to the contour. Otherwise, when the corner point is outside the contour, then the contour closest to the corner point is considered as its subordinate contour.

Fig. 3. Evaluation procedure for affiliation between corner point and contour.
3.2. Trajectory splicing

In some scenarios with high motion velocities, increasing the sampling frequency can effectively improve the measurement accuracy. Shorter time intervals between two frames portray simpler motions of the particles, and provide higher reliability of the optical flow trajectory. If it is necessary to measure particle motion over a longer period of time, the displacement needs to be spliced between frames within the time period.

The trajectory splicing method is shown in Fig. 4. The five points represent the position of a particle a at time t – 2, t – 1, t, t + 1, and t + 2. Sub-track segments, dt – 2, dt – 1, dt, dt + 1, depict the displacement of the particle a at four time intervals. These sub-track segments are highly reliable, because they are matched at smaller time intervals. The velocity vector obtained by connecting at – 2 and at + 2 can be regarded as the correct matching of particles at the interval Δt = (t + 2) – (t – 2), represented by the red line in Fig. 5.

Fig. 4. Reconstruction of particle trajectories over long-time scale.
Fig. 5. Experimental setup. Hopper granular flow width is 15 cm, height is 24 cm, size of the bottom opening is 8 mm, inclination angle of inner wall on both sides of the opening is 30°, and the diameter of the particles is d = 0.8 mm.
4. Experiment and discussion

The section above provides a detailed description of a new PyrLK-based PTV method for complex granular flows. According to the analysis, the measurement accuracy of the method is mainly affected by the selection of the target adjacent region and the number of pyramid layers. The size of the target adjacent region determines how many constraint equations are needed to obtain the optimal solution. The number of equations affects the solution accuracy, which is obtained by the least squares method. The number of Gaussian pyramid layers will determine the measurement accuracy of particles with large displacements. A hopper granular flow device is set up and described in the following section to verify the proposed method.

4.1. Experimental setup

Figure 5 illustrates the experimental environment and measuring device for hopper granular flows. The entire system consists of a hopper device, a light source, a CCD camera, and an image processing device. The hopper granular flow used in this experiment is two-dimensional, with a length of 15 cm, a height of 24 cm, and a thickness of 1 cm. Black glass particles of 0.8-mm size are selected, and the filling process of the particles is distributed evenly.

First, the position of the light source irradiating the hopper is adjusted to ensure that particles receive uniform illumination and form a bright spot. The CCD camera exposure time is set to 5 ms, and images acquired with the camera are processed by a computer. According to the actual image collected by CCD, the number of identifiable particles is about 30000. From the perspective of the number set, 30000 particles can effectively evaluate the performance of the PyrLK-based PTV method for tracking large-scale particle scenarios.

After the particle image is correctly matched, the next step is to select the appropriate parameters. The granular flow velocity at the hopper opening represents the fastest part of the entire granular flow field. Hence, to obtain parameters suitable for the entire hopper particle velocity field, the granular flow at the hopper opening needs to be analyzed. The number of particles in this area is about 1000, and the time interval between two frames is 21 ms. Parameter selection and accuracy values are shown in Table 1. For the hopper granular flow scene, excessive stratification leads to low image quality; a small adjacent region will reduce the number of optical flow constraint equations, and too large adjacent region will contain particles with large velocity contrast. All three cases cause large errors. According to the actual measurement results, this study sets the adjacent region size to 15 × 15 pixels, and the number of layers to two.

Table 1.

Measurement accuracy under different parameters.

.

First, a comparison between the PTV algorithm based on the PyrLK method and the traditional PTV algorithm with matching rules is implemented. These two methods are applied to the hopper granular flow shown in Fig. 5. The measurement area is shown in Fig. 6(a), and the measurement results are shown in Figs. 6(b) and 6(c). The time interval between the two selected matching images is 21 ms.

Fig. 6. (a) Measurement area and velocities of hopper granular flows measured by (b) PyrLK-based PTV and (c) traditional PTV.

The matching error of the traditional PTV is mainly concentrated in the bulk area and outlet, where the velocities exhibit large contrast. According to the principle of PTV, the error of the traditional PTV is caused by corner matching results. To draw a more convincing conclusion, we employ the RANSAC method[29] to detect the mismatching rate of the boundary, bulk, and outlet regions measured by these two methods. The results are listed in Table 2. The mismatching rate of these two methods at the boundary is relatively low, however in the bulk and outlet regions, where the velocities are faster and the contrast is greater, the mismatching rate of the traditional PTV rises sharply, reaching ∼10%, while the error rate of the PyrLK-based PTV method remains within 1%. Moreover, when tracking the motion of 14000 particles, the traditional PTV consumes 0.387 ± 0.001 s, while the PyrLK-based PTV method consumes 0.012 ± 0.001 s, hence the measurement efficiency is increased by an order of magnitude.

Table 2.

Mismatching rate of PyrLK-based PTV and traditional PTV.

.
4.2. Hopper flow results

In summary of the above analysis, the PyrLK-based PTV method is undoubtedly a better choice for measuring large-scale and large-displacement complex particle motion scenarios. Therefore, the following part of this paper involves an experiment to measure and analyze a larger hopper granular flow. In this experiment, motion of the particles lasted for T = 681.673 s, and the initial number of particles in the hopper was ∼30000. Images at four time points were selected for measurement. Figure 7 shows the variation of the velocity field in the granular flow within the hopper. Velocities of particles near the boundary region of the hopper are small, and there is little change at different time. Particles near the outlet area have larger velocities, and considering the bulk area above the outlet, it can be seen that the particle velocities increase from top to bottom.

Fig. 7. Variation of flow velocity field at four moments.

Figure 7 reflects the distribution of the hopper granular flow velocity from a macroscopic perspective. Subsequently, the velocity profiles of typical lateral and longitudinal regions are extracted separately to study the specific motion of the hopper particle.

Figure 8 depicts the velocity distribution of the lateral region at a height of 0.06 m. Velocities of particles at the boundary region of the hopper are very small, and the overall velocities do not increase over time, as the velocity remains around zero. However, as the particles in the hopper become fewer, the velocity of the center position becomes larger. Near the end of the leak, the velocity is almost three times as big as that at the initial time.

Fig. 8. Velocity distribution of lateral region at height of 0.06 m.

Figure 9 depicts the velocity distribution of the vertical position region of the central axis. At any time, the distribution of particle velocities in this region gradually increases from top to bottom, which is evidently caused by the opening of the hopper at the very bottom. In addition, as the particle stacking height is gradually reduced, the particle velocity of the entire central axis increases, and the maximum velocity difference at different time at the same position is close to 3.5 cm/s. However, the particle velocities at the opening are similar, except for the velocity at the initial time.

Fig. 9. Velocity distribution of vertical position region of central axis.
5. Conclusions

In the field of particle motion measurement, because of the existence of numerous complex motion scenes, such as hopper granular flow, the traditional PTV method will lead to a large number of mismatching phenomena, and the measurement efficiency will be reduced accordingly. This study combined PTV technology with an optical flow matching method, improving it by the displacement correction method and trajectory splicing algorithm. The main conclusions of this study are as follows.

(1) Compared with traditional PTV, this improved PTV technology has more accurate matching results and exhibits higher efficiency in the particle measurement field. According to the measurement results of hopper granular flow, the matching error is kept within 1%, and the measurement efficiency is improved by one order of magnitude.

(2) The velocity field distribution of granular flow in a typical hopper is measured. Velocities of particles in the boundary region of the hopper are very low, and as the particles in the hopper become fewer, the velocity of particles at the center position increases. Furthermore, the distribution of particle velocities in this region gradually increases from top to bottom.

Reference
[1] Gong J M Yang H Lin S H Li R Zivkovic V 2018 Powder Technol. 324 76
[2] Schaeper M Damaschke N 2017 Meas. Sci. Technol. 28 055008
[3] Sharp K V Adrian R J 2001 AICHE J. 47 766
[4] Jensen A Pedersen G K 2004 Meas. Sci. Technol. 15 2275
[5] Shi S Ding J Atkinson C Soria J New T H 2018 Exp. Fluids 59 46
[6] Clauser J Knieps M S Büsen M Ding A Schmitz-Rode T Steinseifer U Arens J Cattaneo G 2018 Ann. Biomed. Eng. 46 841
[7] Bolanos-Jimenez R Rossi M Fernandez Rivas D F Kähler C J Marin A 2017 J. Fluid Mech. 820 529
[8] Felix-Felix J R Salinas-Tapia H Bautista-Capetillo C Garcia-Aragon J Burguete J Playan E 2017 Irrig. Sci. 35 515
[9] Ouellette N T Xu H Bodenschatz E 2006 Exp. Fluids 40 301
[10] Maas H G Gruen A Papantoniou D 1993 Exp. Fluids 15 133
[11] Fu S J Biwole P H Mathis C Maissa P 2018 Indoor Built Environ. 27 528
[12] Baek S J Lee S J 1996 Exp. Fluids 22 23
[13] Ferrari G Poletto M 2002 Powder Technol. 123 242
[14] Balevicius R Kacianauskas R Mroz Z Sielamowicz I 2011 Adv. Powder Technol. 22 226
[15] Ma L D Yang G H Zhang S Lin P Tian Y Yang L 2018 Acta Phys. Sin. 67 044501 in Chinese
[16] Zhang S Lin P Yang G H Wan J F Tian Y Yang L 2019 Chin. Phys. B 28 018101
[17] Wang Z W Yang X K Xu Y Yu S Y 2009 Pattern Recognit. Lett. 30 407
[18] Horn B K P Schunck B G 1981 Artif. Intell. 17 185
[19] Ahmine Y Caron G Mouaddib E Chouireb F 2019 Image Vis. Comput. 88 1
[20] Liu Y Xi D G Li Z L Hong Y 2015 J. Hydrol. 529 354
[21] Zhang D J Xie N Liang S Jia J Y 2016 Pattern Recognit. Lett. 76 49
[22] Pinto A M Costa P G Correia M V Matos A C Moreira A P 2017 Robot. Auton. Syst. 87 1
[23] Ohmi K Li H Y 2000 Meas. Sci. Technol. 11 603
[24] Masuda N Ito T Kayama K Kono H Satake S Kunugi T Sato K 2006 Opt. Express 14 587
[25] Qiao Y J Tang Y C Li J S 2013 International Conference on Measurement Information and Control August 16–18, 2013 Harbin, China 1408
[26] Cruz-Santos W Lopez-Garcia L Redondo-Galvan A 2015 Opt. Eng. 54 054102
[27] Barron J L Fleet D J Beauchemin S S 1994 Int. J. Comput. Vis. 12 43
[28] Bouguet J Y 2001 Intel Corporation 5 1
[29] Cao S X Jiang J Zhang G J Yuan Y 2013 Int. J. Remote Sens. 34 2301